library(sf) #for polygons
library(tidyverse) #for data wrangling
library(tmap) #for interactive maps
library(janitor) #for cleaning names
#create list of years
kern_field_years <- c(2016:2019)
#loop through read-in process
for (i in kern_field_years) {
filename <- paste0("kernfields_", i)
wd <- paste0("data/kernfields/kern", i, ".shp")
assign(filename, read_sf(wd))
}
#select region (based off of farmer)
kern_subset <- kernfields_2019 %>%
filter(AGENT == "McMANUS/WILSON,MICHELE/AARIN")
#crop several other years
kern_2017_crop <- st_crop(kernfields_2017, kern_subset)
kern_2018_crop <- st_crop(kernfields_2018, kern_subset)
kern_2019_crop <- st_crop(kernfields_2019, kern_subset)
#create list of cropped fields and tidy
kern_cropped_list <- lapply(mget(sprintf("kern_%d_crop", 2017:2019)),
function(x)
{
clean_names(x) %>%
st_make_valid()
}
)
#bind cropped data together, include new column "source" saying what dataset it came from
sample_bind <- do.call(rbind, lapply(names(kern_cropped_list),
function(x)
cbind(kern_cropped_list[[x]],
source = x)
)) %>%
#turn source column into year column by keeping only numbers
mutate(year = gsub("\\D", "", source))
#merge permit site and year, then select only pmt_site and year
sample_clean <- sample_bind %>%
mutate(pmt_site_year = str_c(year, "_", pmt_site)) %>%
dplyr::select(pmt_site_year) %>%
st_make_valid()
#write to shapefile
st_write(sample_clean, "data/shapefiles_written/sample_clean.shp", append = F)
## Deleting layer `sample_clean' using driver `ESRI Shapefile'
## Writing layer `sample_clean' to data source
## `data/shapefiles_written/sample_clean.shp' using driver `ESRI Shapefile'
## Writing 1410 features with 1 fields and geometry type Unknown (any).
I took the shapefile sample_clean.shp and uploaded it to mapshaper’s online interface (https://mapshaper.org/). Then I entered in the console $ mosaic and downloaded results as a shapefile, as seen in the next code chunk.
This thread could potentially improve workflow https://github.com/mbloch/mapshaper/issues/353
#read in mapshaper output
sample_mosiac <- st_read("data/mapshaper_outputs/sample_clean/sample_clean.shp")
## Reading layer `sample_clean' from data source
## `/Users/irisfoxfoot/Desktop/Kern RA/kernRA_coding_projects/data/mapshaper_outputs/sample_clean/sample_clean.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 817 features and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 6159540 ymin: 2352351 xmax: 6187962 ymax: 2386672
## Projected CRS: NAD83 / California zone 5 (ftUS)
#intersect with field data
sample_join <- st_intersection(st_buffer(sample_mosiac, dist = -1), sample_clean)
#pivot so each column is a year
sample_pivot <- sample_join %>%
separate(pmt_site_year, into = c("year", "pmt_site"), sep = "_") %>%
select(-FID) %>%
pivot_wider(names_from = year, values_from = pmt_site) %>%
st_sf() %>%
st_buffer(dist = 1) %>%
mutate(area_sq_ft = st_area(.))
#unnest listed pmt_site so they are each in own column
sample_pivot_sep <- sample_pivot %>%
unnest_wider(c(`2017`, `2018`, `2019`), names_sep = "_")
#write to csv
write_csv(sample_pivot_sep, "data/shapefiles_written/2017_2019_sample.csv")
#reference what output should be
tmap_mode("view")
tm_shape(sample_bind) +
tm_polygons(id = "pmt_site") +
tm_facets(by = "year", sync = T) +
tm_layout(panel.labels = c("2017", "2018", "2019"))
#test output
tmap_mode("view")
tm_shape(sample_pivot) +
tm_polygons(id = "pmt_site_year") +
tmap_options(check.and.fix = TRUE)